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General conservation equations are derived for 2D dense granular flows within the Boussinesq 
approximation. Assuming that forcing and dissipation equilibrate on average, the relationships be- 
tween local granular temperature, vorticity and stream function are shown to be characterized by 
two functions. Once these two functions are known, the whole velocity and temperature fields can be 
deduced. We checked these predictions on steady surface flows in rotating drum simulated through 
the Non-Smooth Contact Dynamics methods. Also, we briefly discuss some possible ways to derive 
these two functions using statistical mechanics. 
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I. INTRODUCTION 

The intrinsic dissipative nature of the interactions be- 
tween the constituent macroscopic particles sets gran- 
ular media apart from conventional solids, liquids and 
gases [l|. Understanding the rheology of granular sys- 
tems is thus rather difficult. Depending on the flow ve- 
locity, three regimes are usually distinguished: The rapid 
flow - gaseous-like - regime where grains interact through 
binary collisions, is generally described within the frame- 
work of the kinetic theory [2j, bj \M > The slow flow - solid- 
like - regime, where grain inertia is negligible, is most 
commonly described using the tools of soil mechanics and 
plasticity theory 5]. In between these two regimes there 
exists a dense flow - liquid-like - regime where grain in- 
ertia becomes important but contacts between grains are 
still relevant. This last regime has been widely investi- 
gated experimentally, numerically and theoretically (see 
[6( for a review) in various flow configurations. Several 
constitutive laws have been derived by accounting for 
non-local effects @, H, H, El E3 1 by adapting kinetic the- 
ory 0, EH fbil. b y modelling dense flows as partially 
fluidized flows [15j , by considering them as quasi-static 
flows where the mean motion results from transient frac- 
tures modelled as self activated process [H, E3, El, EH 
or more recently by considering them as visco-plastic flu- 
ids [13, HH, ■ To our knowledge, all these approaches 
fail to account for all the features experimentally ob- 
served. In some sense, similar difficulties are encountered 
in the understanding and modelling of turbulent fluids. In 
that case, the challenge is to relate the Reynolds stresses, 
based on small scale fluctuations, to large scale or time 
averaged quantities. A new wa y t o tackle this problem 
was recently suggested [H, [24], |25| , through the consid- 
eration of non-linear steady solutions of the Euler equa- 
tions, thereby disregarding any non-universal effects in- 
duced by (large scale) forcing and (small scale) dissipa- 
tion. When applied to a turbulent von Karman flow, this 



approach leads to the characterization of the steady state 
velocities through two functions, encoding all informa- 
tion about the forcing and the dissipation. In the present 
paper, this method is generalized to the case of inhomo- 
geneous dense granular flows. As a result, one obtains a 
characterization of the steady state through two general 
functions that relate granular temperature, granular vor- 
ticity and granular stream function. In other words, the 
knowledge of these two functions is sufficient to fully re- 
construct the two-dimensional (2D) field of velocity, gran- 
ular temperature, and volume fraction. 



The paper is organized as follows: In Sec. II, the struc- 
ture of 2D steady granular surface flows is derived in 
the rotating drum geometry. Hydrodynamics and state 
equations in granular media are briefly discussed in Sec. 
IIA. Conservation equations are then rewritten assuming 
that volume fraction is nearly constant within the flow 
(Boussinesq approximation) in Sec. IIB, and then in the 
force-free, inviscid limit in Sec. IIC. The general shape of 
the stationary solutions within this limit is given in Sec. 
IID. In particular, it is shown that these states can be 
fully characterized through the knowledge of two func- 
tions F and G. Section III confronts these predictions 
with steady surface flows in rotating drum as obtained 
in Contact Dynamics simulations reported in [26[ that 
were shown to reproduce the experimental features ob- 
served in Refs. [13, [H, H^|. The simulation scheme and 
the description of the simulated systems are briefly re- 
called in Sec. IIIA. Spatial distribution of the averaged 
temperature, volume fraction, vorticity and stream func- 
tion fields are determined within the whole drum, at the 
grain scale (Sec. IIIB). The two characteristic functions 
F and G are then determined from the numerical ex- 
periments (Sec. IIIC) and commented (Sec. HID). In the 
last section of this paper (Sec. IV) some possible ways to 
derive these two functions are briefly discussed. 
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II. THEORETICAL FRAMEWORK: 
CONSERVATION EQUATIONS WITHIN THE 
BOUSSINESQ APPROXIMATION 

A. Granular hydrodynamics 

It is commonly assumed that granular media can be 
described with continuum models. The mass, momentum 
and energy conservation equations then lead to: 

d t v + V • (i/v) = 0, 
d t vw + (v • V) isv = -VP + i/g + F„ sc + F /ore , 
d t isT + V ■ (i/Tv) = -PV ■ v + E msc + E forc . (1) 

In these equations, is(x, z, t) is the field of volume frac- 
tion, v(x, z,t) the velocity field; g is the gravitational 
acceleration; T(x 1 z 1 t) is the field of granular tempera- 
ture defined in term of the RMS part of the velocity 
field, T(x,z,t) = \{{c b (t)-w(x,z,t)) 2 ) b£[x ^ z) where c b (t) 
refers to the instantaneous velocity of the bead b located 
at time t within the elementary volume located in (x, z); 
Ff orc , Ef orc stand for external forcing and F V i SC , E V i SC 
stand for all dissipative processes. This system has to be 
supplemented by an equation of state P — g(y)T and 
rheology, i.e. some constitutive equations describing the 
dissipative processes F„ isc , E visc . 

Contrary to classical liquids, the density and temper- 
ature dependence of transport coefficients play an im- 
portant role in determining the flow density. For dilute 
systems they are usually obtained using kinetic theory 
of granular gases 0, S 0] within the Enskog approxima- 
tion. For dense gases, there is no available systematic the- 
ory allowing their description. They are therefore usually 
prescribed using phenomenological models [1311 or fitted 
using experimental [U H2, [3(| or numerical [3 ll ] data. In 
particular, the equation of state can be written in the 
high-density limit (UHJ: 

P ~ K — T, (2) 

v* — v 

where if is a constant and is* the random closed packed 
limit: is* ~ 0.82 for 2D packing and is* ~ 0.64 for 3D 
packings. At v = j/* this equation therefore predicts a 
zero granular temperature, consistent with the absence 
of motion. 

As for the dissipative terms and forcing, their precise 
shape shall not be needed in the sequel. This is a distin- 
guished feature of our approach. 

B. The Boussinesq approximation 

For simplicity, one focuses on situations where the vol- 
ume fraction is nearly constant close to the random closed 
packed limit v k i/,. In the considered system, this ap- 
proximation is satisfied within 10 percent. Generalization 



to non constant volume fraction is possible, but more in- 
volved. In that limit, the classical Boussinesq approxi- 
mation is implemented by neglecting the fluctuation of 
volume fraction in the continuity equation so that it be- 
comes: 

V-vsiO. (3) 

The other conservation equations may then be simplified 
by defining a reference state with v = 0, T = 0, P = P*, 
v = v*, so that: 

VP* = v*g, (4) 

i.e. an hydrostatic equilibrium in the vertical direction. 
Along with non-zero velocity, we introduce temperature 
and volume fraction deviations with respect to the refer- 
ence state, as: 

v = v* - Sis; T = ST; P = P* + SP. (5) 

The momentum equation can then be written as: 

<V + (v-V)v = --VP + g + F visc + F /orc , 
is 

-i c -i 

~ V5P-— VP* +g VP* 

vise ^forci 

1 Sv 

= VSP g + F„ sc + F/ orc ,(6) 

is* is* 

where the hydrostatic equilibrium has been used to sim- 
plify the last equation. A similar treatment of the tem- 
perature equation leads to: 

P* 

d t 5T+(\r ■ V) ST = E vtsc + Ef orc -V-v « E msc +E forc . 

(7) 

The system of resulting equations can be further trans- 
formed so that it involves only temperature fluctuation 
by using Eq. 

1T = (t-J> (8) 

IS* \ 1 refJ 

where the reference temperature field T re /(r) is given by 
T re f — P*/Kv* 1 so that g = KVT re f, and only the first 
order terms in Sis /is*, 5T/T re f and SP/P* are kept. The 
system of equations of the weakly compressible granular 
medium then takes the shape: 

V • v = 0, 

1 f ST \ 

dtV + (V • V) V = V • SP - ( — ) g + F„ isc + Fforc, 

IS* \ 1 r cfJ 

d t 8T + (v-V)ST= E msc + Eforc- (9) 

Note that the system can also be formulated in a more 
classical Boussinesq-like equation by introducing the 
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variable 6 — ST/T re f and noting that T re f is not a con- 
stant (it varies along the gravity direction), so that: 

V • (v) = 0, 

ftv + (v • V) v = -— V ■ SP -6g + F msc + F forc , 
d t + (v ■ V) 9 + (v • V) logT re/ = E msc + E forc . (10) 

In the sequel, we shall however rather work with the for- 
mulation ©• 



where {ip, <p} = d z ipd x (p — d x ipd z <p is the Jacobian. The 
relation between gravity and T re / was also used to sim- 
plify the buoyancy term. This formulation of the strat- 
ified Navier-Stokes equation has to be supplemented by 
appropriate boundary conditions. Notice that only two 
scalar fields are sufficient to describe the flows under 
consideration: 8T, the granular temperature and q, the 
y-component of the vorticity. 

D. Steady state solutions 



C. 2D case 

We now specialize our granular hydrodynamics to the 
case of 2D medium, such as flow within a thin rotating 
drum of diameter 2R, rotated along the y axis at a con- 
stant angular velocity f2 as investigated in Sec. III. If the 
width of the drum in the y direction is thin with respect 
to the characteristic length scale of (x, z) motions, the 
velocity field can be assumed two-dimensional v(x, z, t). 
In that case, the vorticity is directed along the y axis 
and the forcing is supplied by the boundary conditions. 
In that case, one can recast Eqj9]in cartesian coordinates 
(x, z) as: 



d x v x + d z v z 







(11) 



d t v x + v x d x v x + v z d z v x = — -d x SP - g x ( 

V* \1 ref 

I , rpx 

' vise ' fore ' 

d t v z +v x d x v z +v z d z v z = — -d z SP - g z 

I T?Z _|_ 7712 

' vise r fore 1 

d t ST + v x d x ST + v z d z ST = E msc + Ef orc , 



(12) 



where {v x ,v z ) (resp. (g x ,g z )) denote the components of 
the velocity field (resp. gravity field) in a cartesian ref- 
erential. Note that the third equation expresses the con- 
servation of the granular temperature ST. The two other 
equations for v x and v z involve a pressure field deter- 
mined through incompressibility. However, it can be elim- 
inated by using the stream function ip defined by: 



and 



-d x ip 



The existence of a stream function results from the in- 
compressibility and the 2D nature of the flow. Calling q 
the y-component of the vorticity, one gets: 



d z v x - d x v z = Aip. 



(13) 



where = d x ip + d z ip is the Laplacian. Taking the curl 
of the equation for velocity, the equations (fT2| can be 
recast as: 



d t ST + {ip, ST} 
d t q + q} 



-^visc "I - -Bforc •> (l^) 

K{\ogT ref ,ST} + V x (F 

vise ~\~ F fore) 



Consider now situations where forcing and dissipation 
equilibrate on average: F visc + F forc = E msc + Ef orc = 
, so that one gets steady flows. This situation is the 
one commonly realized in rheometric experiments when 
one imposes a given and tuneable forcing Ff orc , Ef orc to 
get a measured steady velocity profile v. Variations of 
v with respect to fore allows then to find the rheology 
of the studied materials, i.e. the relations between the 
dissipative terms F V i SC , E V i SC and the deformation rate. 

Here, we focus on the left-hand side of Eqs. [14] 
to see the implications of making F V i SC + Ff orc — 
Evisc + Ef orc = on the form taken by the fields ip, q 
and ST. The steady states then obey the averaged equa- 
tions: 



{M = K{\ogT ref ,ST}. 



(15) 



Neglecting fluctuations {ip,ST} « {i/j,6T}, one gets: 

{iP,5T} = 0, (16) 
{V,?} = ^{logT re/ ,<5T}, 

where the over lines over q, T and ip are now omitted for 
sake of simplicity. The first equation is satisfied if 



ST = Fty>), 



(17) 



where F is an arbitrary function. Using the general iden- 
tity 



{f,h(g)} = h'(g){f,g} = {h'(g)f,g}, 



(18) 



Where /, g and h are arbitrary functions, the second 
equation becomes 



{^q + F'^)K\ogT ref } = Q. 



(19) 



Therefore, the general stationary solution of Eqs. (fT4")) is 
of the form 

5T = F{i)) and q + KF'(ip) logT re /= G(ip), (20) 

Where F and G are arbitrary functions. Recalling the 
connection between q and ip, one can therefore charac- 
terize the stationary states through the two function F 
and G as: 

ST = F(ip), 
Aip = q = -KF'{ip)logT ref +G(iP). (21) 
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It should be emphasized that the functions F and G de- 
pend on the forcing and dissipation. Indeed, the com- 
petition between these two effects are responsible for the 
selection of the precise shape for F and G. But once these 
functions are known, one can solve the second equation 
of (f2Tj) to get ip as a function of x and z, and then derive 
from this expression the temperature and velocity profile. 
To close the system of conservation equations, it is then 
sufficient to give the expression for F and G. There are 
probably several ways to prescribe these functions. For 
example, one could use a statistical mechanics approach 
in order to select the "most probable" form of these func- 
tions depending on microscopic constraint s, using meth- 
ods of information theory (see e.g. refs. [23L I24I I32U33I] for 
illustrations in turbulence) . One could also follow the pro- 
cedure used in rheology studies, and try to define these 
functions through "minimal" experimental or numerical 
measurements performed on the considered system. 



III. APPLICATION TO SIMULATED STEADY 
SURFACE FLOWS 

The formalism described in the previous section is now 
applied to the inhomogeneous steady surface flows ob- 
served in rotating drums. 



A. Simulation methodology 



The simulations have 
Smooth Contact Dynamics approach 



been performed using Non- 
" jp. The algo- 
rithms benefit from parallel versions 36j, |37| which show 
their efficiency in the simulation of large systems. The 
scheme has been described in detail elsewhere [26| and 
is briefly recalled below: An immobile drum of diameter 
Dq = 45 cm is half-filled with 7183 rigid disks of density 
po = 2.7 g.cm -3 and diameter uniformly distributed be- 
tween 3 and 3.6 mm. The weak polydispersity introduced 
in the packing prevents 2D ordering effects. The normal 
restitution coefficient between two disks (resp. between 
disks and drum) is set to 0.46 (resp. 0.46) and the friction 
coefficient to 0.4 (resp. 0.95). Once the packing is stabi- 
lized, a constant rotation speed ranging from 2 to 15 rpm 
is imposed to the drum. After one round, a steady contin- 
uous surface flow is reached. One starts then to capture 
400 snapshots equally distributed over a rotation of the 
drum. 

For each bead of each of the 400 frames within a given 
numerical experiment, one records the position r of its 
center of mass and the "instantaneous" velocity c av- 
eraged over the time step St = 6 x 10~ 3 s of the sim- 
ulation. The frame 5ft is defined as the frame rotating 
with the drum that coincides with the reference frame 
5fto = ( e x, e z ) fixed in the laboratory so that e x (resp. e z ) 
is parallel (resp. perpendicular) to the free surface (Fig. 
[T]) . The drum is then divided into elementary square cells 
X(x, z) of size set equal to the mean bead diameter. The 



continuum time-averaged field a(x, z) at a position [x, z) 
is then defined as the average of the corresponding quan- 
tity defined at the grain scale over all the beads in all 
the 400 frames of the sequence whose center of mass is 
within the cell located in (x, z). Figure[T]shows the spatial 
distribution of the x-component of the time-averaged ve- 
locity field v(x, z)&s obtained within this procedure. The 
flowing layer and the static phase are then defined as 
the point where v x is above and below a threshold value 
arbitrary chosen to 0.2. Let us note that all the results 
presented below do not depend on this threhold value. 
The interface between the two phases as defined within 
this procedure is represented as a black line in Fig. [T] 



B. Spatial distribution of the relevant continuous 
fields within the drum 



Let us first determine the granular temperature field 
within the drum. Calling Cj (i) the instantaneous velocity 
of a bead i at a given time t, the fluctuating part of the 
velocity <kij(t) is defined as 5c,(t) = C;(i) — v(x, z) where 
v(x, z) is the mean velocity value on the cell S(x, z) that 
contains the bead i. One can then associate a granular 
temperature Ti(t) = |<Jcf (t) to the considered bead. Fig- 
ure^ shows a typical snapshot of the instantaneous tem- 
perature distribution within the drum as obtained using 
this procedure. Two phases can be clearly distinguished. 
Within the static phase, the temperature is very close 
to zero. Within the flowing layer, the spatial distribution 
of instantaneous local temperature shows large fluctua- 
tions, with hot and cold spots gathered in transient clus- 
ters of various sizes. This structure of hot /cold aggregates 
probably has its origin in the existence of "jammed" ag- 
gregates embedded in the flow, as evidenced in rotating 
drum experiments [381 ] . Since we are primarily interested 
in steady averaged fields in relation with the theoretical 
framework developed in Sec. UU we focus on the spatial 
distribution of the temperature after averaging over the 
400 snapshots of the sequence. The corresponding time- 
averaged temperature field is represented in Fig. ^p. 

Voronoi tessellation is then used to associate an in- 
stantaneous elementary volume as defined in Continuum 
Mechanics to each bead i on each snapshot (see e.g. [26[ 
for related discussion) . Calling Ai the area of the Voronoi 
polyhedra enclosing the grain i, the instantaneous volume 
fraction Vi is defined as Vi — ird 2 /AAi where di denotes 
the diameter of bead i. Typical snapshot of the resulting 
instantaneous map of volume fraction is presented in Fig. 

Apart from a very narrow region - about one bead 
diameter wide - at the free surface and along the drum 
boundary, the volume fraction appears almost constant, 
around (y) ~ 0.82, with apparent random fluctuations 
with standard deviation (v 2 ) — (v) 2 ~ 0.04. However, the 
time-averaged field of volume fraction presented in Fig. 
[5Ji reveals that v decreases slightly within the flowing 
layer, as expected since dilatancy effects should accom- 
pany granular deformation (39l | . 
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To compute the instantaneous vorticity a;, associated 
to each bead i of each snapshot at time t, the following 
procedure is adopted: (i) The Voronoi polygon associated 
with the bead i is dilated homothetically by a factor two, 
so that each edge goes through one of the neighbouring 
beads; (ii) the circulation = J2j c j(t)- s j(t) i s calcu- 
lated around the resulting polygon - each point of a given 
segment Sj is assumed to have a constant velocity Cj (t) 
equal to the one of the embedded bead; (iii) the instanta- 
neous vorticity u^(t) is then defined as qi(t) — Ti(t)/Ai(t) 
where refers to the area of the initial Voronoi poly- 
gon. 

Typical snapshot of the instantaneous vorticity distri- 
bution within the drum as obtained using this procedure 
is presented in Fig. [5^. This distribution is complex. It 
exhibits large fluctuations that self-organize into tran- 
sient network of ID chains. The characterization of this 
transient structure is postponed to future work. Figure 
[2F presents the time averaged vorticity field in the drum 
for Q = 6 rpm. 

The last continuous field of interest in relation with the 
theoretical framework presented in Sec. [TT] is the stream 
function ip(x,z). Its value is set to ip = at the drum 
boundary. The value ip(x, z) is then defined as the flow 
rate going through a line connecting the point M at posi- 
tion r(a;, z) to any point at the drum boundary like e.g. 
point Mo at position ro(:r, — ^/D^/A — x 2 ): 



i/j(x,z)= / v(x, u)v x (x, u)du (22) 

where v(x, z) and v x (x, z) denote the field of volume frac- 
tion and ^-component of velocity computed before. The 
resulting spatial distribution of stream function is shown 
in Fig. O 

C. Determination of the two closure functions 
within Boussinesq approximation 

Let us first determine the value of the parameters K 
and involved in the equation of state given by Eq. 
[21 This determination requires the pressure field P* in 
the reference frame, when v = 0. From Eq. 0] one gets 
P„{x,z) = —zv*cos9 where 9 is the slope of the free 
surface. Equation \E\ can then be rewritten as: 





2 rpm 4 rpm 5 rpm 6 rpm 10 rpm 15 rpm 


K 


0.84 0.86 0.84 0.82 0.83 0.82 
0.4 0.8 1.1 1.1 1 0.8 



TABLE I: Variation of the reference volume fraction v* and 
the parameter K involved in the state Eq. [5] with respect 
to the rotating speed Q. of the drum. They are found to be 
~ 0.83 and K ~ 1 respectively, weakly dependent on the 
rotating speed Q. [iol ]. 

The knowledge of both the field T(x, z) and ip(x, z) al- 
lows to check the first prediction (first equation in system 
HH]) of the approach derived in Sec. [Til Figure [5^ shows 
T(x, z) as a function of ^{x^ z) for = 6 rpm. The points 
clearly gather along a single function. It is worth to em- 
phasize that such result would have been trivial in homo- 
geneous flows such as the observed in Couette geometry 
or inclined plane geometry: in such flows, all the contin- 
uum quantities depend on a single spatial co-ordinate and 
are thus naturally related two by two by single functions. 
On the other hand, the fact that the 2D fields T(x, z) 
and ip(x,z) can be related by a single function in the 
inhomogeneous surface flow considered here, where the 
continuum quantities depend on the both two spatial co- 
ordinates x and z, is highly non trivial and constitutes 
then a rather severe test for the approach derived in Sec. 
HIl The function F (red thick line in Fig. [5^) is then 
defined by averaging the values T falling into logarith- 
mically distributed bins defined along ip. The functions 
F obtained using this procedure for the various rotating 
speeds are represented in Fig. fSb. 

One can now determine the second closure relation 
G(ip). The function F(ip) defined in the previous sec- 
tion (red thick line in Fig. [5^) is first derived numeri- 
cally. The resulting function is then applied in each point 
(x, z) to the field ip(x, z )- Since the reference temperature 
T re f{x, z) = — z cos9/K and the vorticity field q(x, z) are 
also known in each point, one can deduce the value of the 
field q(x, z) + KF'(tp(x, z)) log T re j{x, z) in each point, 
and plot it as a function of ip(x,z) (see Fig. [6^i). The 
points clearly gather along a single function. The func- 
tion G (red thick line in Fig. [5^) is then defined by aver- 
aging the values q(x, z)/KF'(tp(x, z)) \ogT re f{x, z) falling 
into logarithmically distributed bins defined along ip. The 
functions G obtained using this procedure for the various 
rotating speeds f2 are represented in Fig. [SJa. 



v(x, z) = ^ + Kv* , —- (23) 

Z COS V 

The time-averaged local volume fraction v{x, z) is plot- 
ted as a function of the ratio T(x, z)/z cos 9 in Fig. [U The 
values of both ^» and K can then be deduced (see TabQ}. 
They are found to be ~ 0.83 and K ~ 1 respectively, 
weakly dependent on the rotating speed |4(J. The ref- 
erence temperature field T re f(x, z) = —zcos9/K is then 
known. 



D. Discussion of the results 

Our determination of the two closure functions calls 
for some comments. A first noticeable feature is that 
the closure function extends smoothly, without any no- 
ticeable transition, from the static to the flowing region. 
This is quite remarkable, since both phases are character- 
ized by different dynamical properties, and since our hy- 
drodynamic description presumably applies best within 
the flowing region. The main difference between the two 
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phases is in the scattering of the data along the fit: It is 
larger in the flowing region than in the static region. This 
may be traced to correlated fluctuations that have been 
neglected in our approach, see after Eq. 15, and that are 
larger in the flowing region. It would be interesting to 
see if a longer time averaging leads to a reduction of this 
scattering. 

An interesting comparison can also be made with re- 
spect to a real fluid system, where a similar approach can 
be used and where dissipation is made through ordinary 
viscosity In that case, it has been shown in [25| that the 
determination of the closure function is valid only in the 
bulk flow region. Outside this region, the data scatters 
randomly, without forming any specific shape. A possi- 
ble explanation was that outside the bulk, i.e. closer to 
the boundaries and the flow forcing devices, viscous and 
forcing processes become important and do not balance 
locally on average as assumed here. The reason why it 
works so well in the granular case, without any need of 
selecting any flow region, may lie in the local character 
of the dissipative processes that precludes any long-range 
correlation between forcing and dissipation. 

IV. CONCLUDING DISCUSSION 

In this paper, the stationary states of a rotating dense 
inhomogeneous granular flow have been shown to be 
characterized by two functions F and G linking the gran- 
ular temperature, vorticity and stream functions. When 
varying the rotation speed of the apparatus, both F and 
G vary, so that they can be seen as representative of the 
state of the system. Conversely, the knowledge of F and 
G is sufficient to fully reconstruct the 2D fields for tem- 
perature and velocity. Let us first consider the flowing 
phase. From Fig. \5\ and Fig. [6l one sees that, in that 
phase, F is asymptotically linear F ~ aip, while G is ap- 
proximately constant, G ~ b. Inserting these shapes in 
cq. (I2ip leads to: 

ST = atp, 

Aip = -aK log(- z cos 6 /K)+b (24) 

Integrating the second equation with respect to z, one 
finds: 

z 2 3 aK 
tp = (b- aKlog(cos8/K)) — + -aKz —z 2 \og(-z) 



so that the temperature profile is quadratic, with loga- 
rithmic correction and the velocity profile is linear, with 
logarithmic correction. This is indeed the behaviour ob- 
served in this type of flow in the flowing phase [f| [ll], [26|, 

HI. 

In the static phase, F appears quadratic in tp , F ~ cip 2 
and G is linear G ~ dip. Therefore, eq. (f2"Tj) becomes: 

5T = cip 2 , 

Aip = (d-2cKlogT ref )ip (26) 
The solution for ip is in this case 

tp = tp exp(h(z)), 
h' 2 (z) + h"(z) = d- 2cK\ogT ref > 0, (27) 

so that both the velocity profile and the temperature 
profiles are exponential, with algebraic corrections. This 
is indeed the behaviour observed in the static phase of 
the rotating drum [28l. lill. [42l| . 

These simple asymptotics show the importance of un- 
derstanding the shape of F and G as a function of the 
forcing and dissipation. From an experimental or numer- 
ical point of view, one may try and find empirical laws 
from variation of the control parameters like rotation 
speed, size of the beads, friction coefficient, etc. From 
a theoretical point of view, one could try to link these 
functions to the macroscopic flow geometry and the mi- 
croscopic properties of the flow or rheology through ki- 
netic equations. In a next paper, we explore an alternate 
strategy, based upon the statistical mechanics Jaynes' in- 
terpretation in term of information theory. This will lead 
to Gibbs distributions providing a direct link between 
the function F and G and the fluctuations of physical 
quantities. In this respect, the present approach provides 
a useful insight for dense granular flow. 
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FIG. 1: Time averaged velocity field in the simulated 2D rotating drum. In this peculiar numerical experiment, the rotation 
speed was S7 = 6 rpm. The black line shows the interface between the flowing layer and the "static" packing. Velocities are 
expressed in m/s (left-hand colorbar) or non-dimensionalized by \/~(gd) (right-hand colorbar) where g refers to the gravity 
constant and d to the mean diameter of the beads (see text for details). 
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FIG. 2: Spatial dististribution of the various continuous fields measured within the drum. In this peculiar numerical experiment, 
the rotation speed was fi = 6 rpm. Left: Typical snapshot of the instantaneous spatial distribution of granular temperature 
(a), volume fraction (c) and vorticity (e) within the rotating drum. Right: Time averaged field of temperature (b), volume 
fraction (d) and vorticity (f). The average was taken over the 400 snapshots of the sequence for $1 = 6 rpm. Temperatures are 
expressed in m 2 /s 2 (left-hand colorbar) or non-dimensionalized by gd (right-hand colorbar). Vorticities are expressed in s _1 
(left-hand colorbar) or non-dimensionalized by y/gjd (right-hand colorbar) where g refers to the gravity constant and d to the 
mean diameter of the beads (see text for details) . 
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FIG. 3: Spatial distribution of the stream function ip in the drum for £7 = 6 rpm. The stream function is expressed in m 2 /s 
(left-hand colorbar) or non-dimensionalized by d^J\gd) (right-hand colorbar) where g refers to the gravity constant and d to 
the mean diameter of the beads (see text for details). 
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FIG. 4: Volume fraction v(x, z) as a function of the ratio T(x, z)/zcos9 for = 6 rpm. For this rotation speed, the mean slope 
of the free surface was measured to be ~ 19.7° [26J. The straight line is a fit given by Eq. [53] with K ~ 1.1 and v* ~ 0.82. 
The temperature is non-dimensionalized hygd where g refers to the gravity constant and d to the mean diameter of the beads 
(see text for details). 
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FIG. 5: (a): Variation of the local temperature T(x, z) as a function of the local stream function ip(x, z) for Q, = 6 rpm. Each 
dark dot of the cloud corresponds to an elementary square cell ~S(x, z) of size equal to the mean bead diameter. Green dots 
(resp. blue dots) correspond to points that belong to the flowing layer (resp. to the static phase). The red thick line shows the 
function (T) = F(ip) where (T) is defined as the average of the values T of the cells z) whose if>(x, z) fall into logaritmically 
distributed bins, (b) Variation of F(ip) as a function of Q. The temperature and stream function are non-dimensionalized by 
gd and d^fgd respectively, where g refers to the gravity constant and d to the mean diameter of the beads (see text for details). 
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FIG. 6: (a): Variation of the local field q{x,z) + KF'(4>(x, z)) log T re f(x, z) as a function of the local stream function ip(x,z) 
for Q = 6 rpm. Each dark dot of the cloud corresponds to an elementary square cell z) of size equal to the mean bead 
diameter. Green dots (resp. blue dots) correspond to points that belong to the flowing layer (resp. to the static phase). The red 
thick line shows the function {q + KF'(tp) logT re /) = G{ip) where () is defined as the average on the cells E(a;, z) with values 
tp(x, z) which fall into logaritmically distributed bins, (b) Variation of G(ip) as a function of fi. The temperature, vorticity and 
stream function are non-dimensionalized by gd, \J g/d and d^fgd respectively, where g refers to the gravity constant and d to 
the mean diameter of the beads (see text for details). 



